clear all; close all; clc;
location = strcat(pwd,'\');
varin = strcat(location,'data\');
figout = strcat(location,'figures\');
path(path,strcat(location,'function_calls'))

% Load value function across parameter space
ann_alpha = 3:0.3:9; % true alpha [3 to 9 by 1s]
ann_sigma = 11:0.8:19; % signal precision [aum+returns]

n_alpha = cell(size(ann_alpha,2),1); % depository for intensities
vf_alpha = cell(size(ann_alpha,2),1); % depository for value function
pi_alpha = cell(size(ann_alpha,2),1); % depository for inv value
for ali=1:size(ann_alpha,2)
    alpha = ann_alpha(ali);
    
    ind = (size(ann_sigma,2)+1)/2;
    sigma = ann_sigma(ind); % mid sigma
    
    load(strcat(varin,'vf','_',num2str(alpha*10),'_',num2str(sigma*10),'.mat'));

    % Define absolute intensity measure, sigma^2*n/(1+n)
    nName = strcat('_',num2str(alpha*10),'_',num2str(sigma*10));    
    n_alpha{ali} = eval(strcat('var',nName))*eval(strcat('inten',nName))./...
        (1+eval(strcat('inten',nName)));
    vf_alpha{ali} = eval(strcat('v',nName));
    pi_alpha{ali} = eval(strcat('pi',nName));
end

n_sigma = cell(size(ann_sigma,2),1);
vf_sigma = cell(size(ann_sigma,2),1);
pi_sigma = cell(size(ann_sigma,2),1);
for sigj=1:size(ann_sigma,2) 
    ind = (size(ann_alpha,2)+1)/2; 
    alpha = ann_alpha(ind); % mid alpha
    
    sigma = ann_sigma(sigj); 
    
    load(strcat(varin,'vf','_',num2str(alpha*10),'_',num2str(sigma*10),'.mat'));
    
    % Define absolute intensity measure, sigma^2*n/(1+n)
    nName = strcat('_',num2str(alpha*10),'_',num2str(sigma*10));    
    n_sigma{sigj} = eval(strcat('var',nName))*eval(strcat('inten',nName))./...
        (1+eval(strcat('inten',nName)));    
    vf_sigma{sigj} = eval(strcat('v',nName));    
    pi_sigma{sigj} = eval(strcat('pi',nName));    
end

%% Contour Plots, alpha = 6%/yr & sigma = 15%/yr
fig2A = figure;
nMdls = 2;
Axes = zeros(nMdls);
Axes(1) = axes;

% Investment Decision
ind = (v_60_150==pi_60_150 & p>0.1 & q<0.9 & p-q>0.25);
boundaries = [1,1];
[~,h] = contour(q,p,ind,boundaries,'k-.');

h.LineWidth = 4;
set(gca,'visible','on')
xlabel('Market Assess. (P(High | Public))')
ylabel('Allocator''s Assess. (P(High | Public + Private))')
ytickformat('%.1f')
xtickformat('%.1f')

ax = gca;
ax.FontSize = 14;
hold on;

ind = (p==q);
temp_x = p(ind);
temp_y = q(ind);
plot(temp_x,temp_y,'k:')

text(0.05, 0.87,'Investment Region','FontSize',16)
text(0.3,0.5,'Due-diligence Region','FontSize',16)

set(fig2A,'PaperSize',fliplr(get(fig2A,'PaperSize'))) 
print(fig2A,'-fillpage','-dpdf','-r500',strcat(figout,'Figure3A'))

% Intensity Decision
fig2B = figure;
nMdls = 2;
Axes = zeros(nMdls);
Axes(1) = axes;

% Investment Decision
ind = (v_60_150==pi_60_150 & p>0.1 & q<0.9 & p-q>0.25);
boundaries = [1,1];
[~,h] = contour(q,p,ind,boundaries,'k-.');
h.LineWidth = 4;
set(gca,'visible','on')
xlabel('Market Assess. (P(High | Public))')
ylabel('Allocator''s Assess. (P(High | Public + Private))')
ytickformat('%.1f')
xtickformat('%.1f')

ax = gca;
ax.FontSize = 14;
hold on;

ind = (p==q);
temp_x = p(ind);
temp_y = q(ind);
plot(temp_x,temp_y,'k:')
hold on;

n_60_150 = inten_60_150; % define absolute intensity
ind = (p>0.98);
n_60_150(ind) = 0;

ns = n_60_150(:); 
nmax = max(ns(:));
nmin = min(ns(:));

inc = 20; % increment contour lines
ninc = (nmax - nmin)/inc; 
nlvl = nmin:ninc:nmax;
nlvl = nlvl(12:end);

% Fix gray scale across panel
map = gray;
map = flipud(map);
map = map(1:floor(size(map,1)/inc):size(map,1),:);

Axes(2) = axes;
colormap(map);

contour(q,p,n_60_150,nlvl);

% Hide the top axes
set(Axes(2), 'visible', 'off');
set(Axes(2), 'XTick', []);
set(Axes(2), 'YTick', []);

set(fig2B,'PaperSize',fliplr(get(fig2B,'PaperSize'))) 
print(fig2B,'-fillpage','-dpdf','-r500',strcat(figout,'Figure3B'))

linkaxes(Axes)
hold on;

% Probable priors
rectangle('position',[0.01 0.05 0.30 0.70])
rectangle('position',[0.50 0.75 0.45 0.24])
hold on;

% Path of highest gradient
quiver(0.03,0.07,0.27,0.61,0,'--k','MaxHeadSize',0.5) % arrow from (0.03,0.07) to (0.30, 0.68)
quiver(0.93,0.97,-.38,-.12,0,'--k','MaxHeadSize',0.5) % arrow from (0.03,0.07) to (0.30, 0.68)
hold on;

% Hide the top axes
set(Axes(2), 'visible', 'off');
set(Axes(2), 'XTick', []);
set(Axes(2), 'YTick', []);

set(fig2B,'PaperSize',get(fig2B,'PaperSize')) 
print(fig2B,'-fillpage','-dpdf','-r500',strcat(figout,'Figure3D'))

%% Intensity vs Spread
fig3A = figure;

inten_func = griddedInterpolant(q',p',inten_60_150');

% Panel A cross section
q_quiver_young = (0.07:0.001:0.30)'; % see above, assume linear relationship along line of interest
q_quiver_old = (0.85:-0.001:0.55)';

% assume linear --> y = mx + b, solve for m and b
b = (0.68-(0.30*0.07/0.03))/(1-0.27/0.03);
m = (0.68-b)/0.30;
p_quiver_young = m*q_quiver_young + b;

b = (0.97-(0.93*0.85/0.55))/(1-0.93/0.55);
m = (0.97-b)/0.93;
p_quiver_old = m*q_quiver_old + b;

lin_young = inten_func(q_quiver_young,p_quiver_young);
s_young = p_quiver_young-q_quiver_young;

lin_old = inten_func(q_quiver_old,p_quiver_old);
s_old = p_quiver_old-q_quiver_old;

% Plot
lin_young = (lin_young-min(lin_young))/max(lin_young-min(lin_young)); % normalize, fraction of max
lin_old = (lin_old-min(lin_old))/max(lin_old-min(lin_old));

plot(s_young,lin_young,'k','LineWidth',4); hold on;
plot(s_old,lin_old,'k--','LineWidth',4)

ylim([0.0 1.05]); xlim([0.05 0.40]);
xlabel('Allocator minus Market Assessment')
ylabel('Max Intensity')
yticks([0.0 0.15 0.30 0.45 0.60 0.75 0.90 1.05])
xtickformat('%.2f')
ytickformat('%.2f')

text(0.15, 0.92,'Established','FontSize',16)
text(0.26,0.45,'Young','FontSize',16)

ax = gca;
ax.FontSize = 16;
set(fig3A,'PaperSize',fliplr(get(fig3A,'PaperSize'))) 
print(fig3A,'-fillpage','-dpdf','-r500',strcat(figout,'Figure4A'))

%% Max Intensity vs Alpha
fig3C = figure;

% Compute max intensity for each vf (ignore activity at very high p levels)
inten_vs_alpha = nan(size(n_alpha));

for i=1:size(n_alpha,1)
    inten = n_alpha{i};

    ind = (p>0.98);
    inten(ind) = 0;

    for j = 1:101
        inten(:,j) = smooth(inten(:,j));
    end

    inten = inten';
    
    for j = 1:101
        inten(:,j) = smooth(inten(:,j));
    end
    
    inten = inten';

    inten_vs_alpha(i) = max(max(inten));
end

% Plot
max_n = max(inten_vs_alpha);
inten_vs_alpha = inten_vs_alpha/max_n;

ann_csigma = (1/3)*ann_alpha.^2; % dispersion assuming uniform density

plot(ann_csigma,smooth(inten_vs_alpha),'k','LineWidth',4)

ylim([0.8 1.05]); xlim([0 30]);
xlabel('Fund Quality Dispersion')
ylabel('Max Intensity')
yticks([0.8 0.85 0.9 0.95 1.0 1.05])
xtickformat('%.1f')
ytickformat('%.2f')

ax = gca;
ax.FontSize = 16;
set(fig3C,'PaperSize',fliplr(get(fig3C,'PaperSize'))) 
print(fig3C,'-fillpage','-dpdf','-r500',strcat(figout,'Figure4C'))

%% Max Intensity vs Sigma
fig3B = figure;

% Compute max intensity for each vf
inten_vs_sigma = nan(size(n_sigma));
for i=1:size(n_sigma,1)
    inten = n_sigma{i};
    vf = vf_sigma{i};
    pi = pi_sigma{i};
    
    ind = (vf==pi & p>0.1 & q<0.9 & p-q>0.25);
    boundaries = [1,1];
    [c,h] = contour(q,p,ind,boundaries,'k-.');

    lev = c(1,1);
    cnt = c(2,1);

    inten_func = griddedInterpolant(q',p',inten');
    inten_inv_line = inten_func(0.005+c(1,3:cnt),c(2,3:cnt)-0.005); % intensity along investment line, right before investment

    inten_vs_sigma(i) = max(inten_inv_line);
end

% Plot
max_n = max(inten_vs_sigma);
inten_vs_sigma = inten_vs_sigma/max_n;

ann_info = (1./(ann_sigma/100).^2); % precision or "informativeness"

plot(ann_info,smooth(inten_vs_sigma),'k','LineWidth',4)

ylim([0.55 1.05]); xlim([20 90]);
xlabel('Signal Informativeness')
ylabel('Max Intensity')
yticks([0.55 0.65 0.75 0.85 0.95 1.05])
xtickformat('%.1f')
ytickformat('%.2f')

ax = gca;
ax.FontSize = 16;
set(fig3B,'PaperSize',fliplr(get(fig3B,'PaperSize'))) 
print(fig3B,'-fillpage','-dpdf','-r500',strcat(figout,'Figure4B'))